#include "Element_plane62.h"

#include "Node_Base.h"

void Element_plane62::hammer_integral2D()
{
    ci.resize(4);
    si.resize(4);
    wi.resize(4);

    ci[0] = 1.0 / 3;
    ci[1] = 0.6;
    ci[2] = 0.2;
    ci[3] = 0.2;

    si[0] = 1.0 / 3;
    si[1] = 0.2;
    si[2] = 0.6;
    si[3] = 0.2;

    wi[0] = -27.0 / 48 * 0.5;
    wi[1] = 25.0 / 48 * 0.5;
    wi[2] = 25.0 / 48 * 0.5;
    wi[3] = 25.0 / 48 * 0.5;
}

void Element_plane62::exhammer()
{
    ci.resize(6);
    si.resize(6);
    zi.resize(6);

    ci[0] = 2;
    ci[1] = -0.5;
    ci[2] = -0.5;
    ci[3] = 0;
    ci[4] = 0.75;
    ci[5] = 0.75;

    si[0] = -0.5;
    si[1] = 2;
    si[2] = -0.5;
    si[3] = 0.75;
    si[4] = 0;
    si[5] = 0.75;

    zi[0] = -0.5;
    zi[1] = -0.5;
    zi[2] = 2;
    zi[3] = 0.75;
    zi[4] = 0.75;
    zi[5] = 0;
}

void Element_plane62::Nodeff()
{
    ci.resize(6);
    si.resize(6);

    ci[0] = 1;
    ci[1] = 0;
    ci[2] = 0;
    ci[3] = 0;
    ci[4] = 0.5;
    ci[5] = 0.5;

    si[0] = 0;
    si[1] = 1;
    si[2] = 0;
    si[3] = 0.5;
    si[4] = 0;
    si[5] = 0.5;
}

void Element_plane62::get_N(double c, double s, Eigen::VectorXd& N)
{
    N.resize(6);

    const double z = 1 - c - s;

    N(0) = c * (2 * c - 1);
    N(1) = s * (2 * s - 1);
    N(2) = z * (2 * z - 1);
    N(3) = 4 * s * z;
    N(4) = 4 * c * z;
    N(5) = 4 * c * s;
}

void Element_plane62::get_Nex(double c, double s, double z, Eigen::VectorXd& N)
{
    //Eigen::VectorXd Nten;
    //Nten.resize(10);

    //Nten(0) = 0.5 * c * (3 * c - 1) * (3 * c - 2);
    //Nten(1) = 0.5 * s * (3 * s - 1) * (3 * s - 2);
    //Nten(2) = 0.5 * z * (3 * z - 1) * (3 * z - 2);
    //Nten(3) = 4.5 * c * s * (3 * c - 1);
    //Nten(4) = 4.5 * c * s * (3 * s - 1);
    //Nten(5) = 4.5 * s * z * (3 * s - 1);
    //Nten(6) = 4.5 * s * z * (3 * z - 1);
    //Nten(7) = 4.5 * z * c * (3 * z - 1);
    //Nten(8) = 4.5 * z * c * (3 * c - 1);
    //Nten(9) = 27.0 * c * s * z;

    //Eigen::MatrixXd M;
    //M.resize(10, 4);
    //M.setZero();
    //M(0, 0) = 1.0;
    //M(1, 1) = 1.0;
    //M(2, 2) = 1.0;
    //M(3, 0) = 2.0 / 3;
    //M(3, 1) = 1.0 / 3;
    //M(4, 0) = 1.0 / 3;
    //M(4, 1) = 2.0 / 3;
    //M(5, 1) = 2.0 / 3;
    //M(5, 2) = 1.0 / 3; 
    //M(6, 1) = 1.0 / 3;
    //M(6, 2) = 2.0 / 3;
    //M(7, 2) = 2.0 / 3;
    //M(7, 0) = 1.0 / 3;
    //M(8, 2) = 1.0 / 3;
    //M(8, 0) = 2.0 / 3;
    //M(9, 3) = 1.0;
    //
    //N = Nten.transpose() * M;

    N.resize(3);
    N(0) = c;
    N(1) = s;
    N(2) = z;
}

void Element_plane62::get_dNcs(double c, double s, Eigen::VectorXd& dNc, Eigen::VectorXd& dNs)
{
    dNc.resize(6);
    dNs.resize(6);

    const double z = 1 - c - s;

    dNc(0) = 4 * c - 1;
    dNc(1) = 0;
    dNc(2) = 1 - 4 * z;
    dNc(3) = -4 * s;
    dNc(4) = 4 * (z - c);
    dNc(5) = 4 * s;

    dNs(0) = 0;
    dNs(1) = 4 * s - 1;
    dNs(2) = 1 - 4 * z;
    dNs(3) = 4 * (z - s);
    dNs(4) = -4 * c;
    dNs(5) = 4 * c;
}

double Element_plane62::get_csyz(double c, double s, Eigen::Matrix2d& cs_yz)
{
    Eigen::VectorXd dNc;
    Eigen::VectorXd dNs;
    get_dNcs(c, s, dNc, dNs);
    double y_c = 0, y_s = 0, z_c = 0, z_s = 0;
    for (int i = 0; i < 6; i++)
    {
        auto pNode = m_pNode[i].lock();
        y_c += dNc(i) * pNode->m_Coord[0];
        y_s += dNs(i) * pNode->m_Coord[0];
        z_c += dNc(i) * pNode->m_Coord[1];
        z_s += dNs(i) * pNode->m_Coord[1];
    }
    Eigen::Matrix2d yz_cs;
    yz_cs << y_c, y_s, z_c, z_s;
    double det = yz_cs.determinant();
    cs_yz = yz_cs.inverse();

    return det;
}

double Element_plane62::get_dNyz(double c, double s, Eigen::VectorXd& dNy, Eigen::VectorXd& dNz)
{
    Eigen::Matrix2d cs_yz;
    double det = get_csyz(c, s, cs_yz);

    Eigen::VectorXd dNc;
    Eigen::VectorXd dNs;
    get_dNcs(c, s, dNc, dNs);
    double& c_y = cs_yz(0, 0);
    double& c_z = cs_yz(0, 1);
    double& s_y = cs_yz(1, 0);
    double& s_z = cs_yz(1, 1);

    dNy.resize(6);
    dNz.resize(6);
    for (int i = 0; i < 6; i++)
    {
        dNy(i) = dNc(i) * c_y + dNs(i) * s_y;
        dNz(i) = dNc(i) * c_z + dNs(i) * s_z;
    }

    return det;
}

double Element_plane62::get_y(Eigen::VectorXd& N)
{
    double y;
    Eigen::VectorXd yi;
    yi.resize(6);
    for (int i = 0; i < 6; i++)
    {
        yi(i) = m_pNode[i].lock()->m_Coord[0];
    }

    y = yi.dot(N);

    return y;
}

double Element_plane62::get_z(Eigen::VectorXd& N)
{
    double z;
    Eigen::VectorXd zi;
    zi.resize(6);
    for (int i = 0; i < 6; i++)
    {
        zi(i) = m_pNode[i].lock()->m_Coord[1];
    }

    z = zi.dot(N);

    return z;
}
